!======================================================================
  subroutine particle_move
! Move particles every time step
! G. Ito 3/07
!======================================================================

include 'precision.inc'
include 'params.inc'
include 'arrays.inc'

integer iorder(np)
!real*8 xpord(np)
!--------------------------------------------------------------------------
!Interpolate velocities to particles and move them
!--------------------------------------------------------------------------
dxp=rxbo/dble(np)

do 100 i=1,np
  vxp=999.
  vzp=999.
  if (xp(i).lt.cord(1,1,1)) then		!particle outside left side of nodes
    vxp=vel(1,1,1)
    vzp=vel(1,1,2)
  elseif (xp(i).gt.cord(1,nx,1)) then		!particle outside right side of nodes
    vxp=vel(1,nx,1)
    vzp=vel(1,nx,2)
  else						!particle within bounds of nodes
    do 110 ii=1,nx-1
      xl=cord(1,ii,1)
      xr=cord(1,ii+1,1)
      if (xp(i).ge.xl.and.xp(i).le.xr) then
        vxl=vel(1,ii,1)
        vxr=vel(1,ii+1,1)
        vzl=vel(1,ii,2)
        vzr=vel(1,ii+1,2)
        vxp=vxl + ((vxr-vxl)/(xr-xl)) * (xp(i)-xl)
        vzp=vzl + ((vzr-vzl)/(xr-xl)) * (xp(i)-xl)
        exit
      endif
    110 continue
  endif
  
  xp(i)=xp(i) + vxp*dt
  zp(i)=zp(i) + vzp*dt
  
  if (vxp.eq.999.or.vzp.eq.999.0) then
    write(*,*) 'Problems in particle_move, setting particle velocity'
    stop
  endif
!  write(*,*) xp(i), zp(i)
100  continue


!--------------------------------------------------------------------------
!Re-seed if particles outside of model domain
!-------------------------------------------------------------------------- 
do 300 ip=1,np
  if ( (xp(ip).gt.(x0+rxbo+dxp)) .or. (xp(ip).lt.(x0-dxp))) then
    call sortrx(np,xp,iorder)		!identify order of increasing xp
    dxpmax=-99.0
    zp(ip)=9999.0
    do 200 ipp=1,np-1			!Find midpoint between to most distal particles
      ip1=iorder(ipp)
      ip2=iorder(ipp+1)
      dxp=dabs(xp(ip2)-xp(ip1))
      if (dxp .gt. dxpmax.and.(ip1.ne.ip).and.(ip2.ne.ip)) then
        dxpmax=dxp
        ipmax1=ip1
	ipmax2=ip2
      endif
!      write(*,'(4i4,4f10.2)') ip, ipp,ip1,ip2, xp(ip)/1e3, xp(ip1)/1e3, xp(ip2)/1e3,dxpmax
    200 continue  

    xp(ip)=(xp(ipmax2)+xp(ipmax1))/2.0  
    
    do 250 ii=1,nx-1 			!Now find z-coord of new particle position
      xl=cord(1,ii,1)
      xr=cord(1,ii+1,1)
      if (xp(ip).ge.xl.and.xp(ip).le.xr) then
        zl=cord(1,ii,2)
        zr=cord(1,ii+1,2)
        zp(ip)=zl + ((zr-zl)/(xr-xl)) * (xp(ip)-xl)
        exit
      endif
    250 continue
    if (zp(ip).eq.9999.0) then
      write(*,*) 'Problems re-seeding in particle_move, ip, xp, zp=', ip, xp(ip), zp(ip)
      stop
    endif
  endif
300 continue

return
end

